#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "papi.h" /* This needs to be included every time you use PAPI */

#define ERR -1
#define IJK 1
#define JIK 2
#define JKI 3
#define KJI 4
#define IKJ 5
#define KIJ 6
#define NUM_EVENTS 1
#define ERROR_RETURN(retval) { fprintf(stderr, "Error %d %s:line %d: \n", retval,__FILE__,__LINE__);    exit(retval); }

int atoloop_order(char* arg);
double** alloc_matrix(int size);
void cache_reset();
void free_matrix(double** M, int size);
void matrix_ones_init(double** M, int size);
void matrix_rand_init(double** M, int size);
void matrix_zero_init(double** M, int size);
void print_matrix(double** M, int size);

int main(int argc, char** argv)
{
    int retval;

    /* Must be initialized to PAPI_NULL before calling PAPI_create_event */
    int EventSet = PAPI_NULL;

    char errstring[PAPI_MAX_STR_LEN];
    long long values[NUM_EVENTS];

    if (argc != 3)
    {
        fprintf(stdout, "Usage: %s [size] [loop order: ijk|jik|jki|kji|ikj|kij]\n", argv[0]);

		return 1;
    }

    int size = atoi(argv[1]);
	int loop_order = atoloop_order(argv[2]);

	if (loop_order == ERR)
	{
		fprintf(stdout, "Invalid loop order.\n\n");
		fprintf(stdout, "Usage: %s [size] [loop order: ijk|jik|jki|kji|ikj|kij]\n", argv[0]);
		
		return 1;
	}

    /* Allocate matrices spaces */
    double** A = alloc_matrix(size);
    double** B = alloc_matrix(size);
    double** C = alloc_matrix(size);

    /* Initialize matrices */
    matrix_ones_init(A, size);
    matrix_ones_init(B, size);
    matrix_zero_init(C, size);

    int i, j, k;

    cache_reset();
	
	/****************************************************************************
     * This part initializes the library and compares the version number of the *
     * header file, to the version of the library, if these don't match then it *
     * is likely that PAPI won't work correctly. If there is an error, retval   *
     * keeps track of the version number.                                       *
     ****************************************************************************/
    if((retval = PAPI_library_init(PAPI_VER_CURRENT)) != PAPI_VER_CURRENT )
    {
        fprintf(stderr, "Error: %s\n", errstring);
        exit(1);
    }

    /* Creating event set */
    if ((retval=PAPI_create_eventset(&EventSet)) != PAPI_OK)
        ERROR_RETURN(retval);
    
    /* Add the array of events PAPI_TOT_INS and PAPI_TOT_CYC to the eventset */
    if ((retval=PAPI_add_event(EventSet, PAPI_L2_DCA)) != PAPI_OK)
        ERROR_RETURN(retval);
    
    /* Start counting */
    if ( (retval=PAPI_start(EventSet)) != PAPI_OK)
        ERROR_RETURN(retval);
        
	/* Matrix-Matrix Multiply Calculation */
	// Watch out when permuting the loops by renaming variables.
    for (i=0; i<size; i++)
    {
        for (j=0; j<size; j++) 
        {
            for (k=0; k<size; k++)
            {
				if (loop_order == IJK)
				{
					C[i][j] = C[i][j] + A[i][k] * B[k][j];
				}
				//       loop_order == IJK
				else if (loop_order == JIK)
				{
					C[j][i] = C[j][i] + A[j][k] * B[k][i];
				}// C[i][j] = C[i][j] + A[i][k] * B[k][j];
				//       loop_order == IJK
				else if (loop_order == JKI)
				{
					C[k][i] = C[k][i] + A[k][j] * B[j][i];
				}// C[i][j] = C[i][j] + A[i][k] * B[k][j];
				//       loop_order == IJK
				else if (loop_order == KJI)
				{
					C[k][j] = C[k][j] + A[k][i] * B[i][j];
				}// C[i][j] = C[i][j] + A[i][k] * B[k][j];
				//       loop_order == IJK
				else if (loop_order == IKJ)
				{
					C[i][k] = C[i][k] + A[i][j] * B[j][k];
				}// C[i][j] = C[i][j] + A[i][k] * B[k][j];
				//       loop_order == IJK
				else if (loop_order == KIJ)
				{
					C[j][k] = C[j][k] + A[j][i] * B[i][k];
				}// C[i][j] = C[i][j] + A[i][k] * B[k][j];
            }
        }
    }

    /* Stop counting, this reads from the counter as well as stop it. */
    if ( (retval=PAPI_stop(EventSet,values)) != PAPI_OK)
        ERROR_RETURN(retval);

    double tot_access = 4.0 * (double)size * (double)size * (double)size;
	double cache_misses = (double)values[0];
    double miss_ratio = cache_misses / tot_access;
    printf("%s\t%d\t%f\t%f\t%f\n", argv[2], size, miss_ratio, cache_misses, tot_access);
    
    if ( (retval=PAPI_remove_event(EventSet,PAPI_L2_DCA))!=PAPI_OK)
        ERROR_RETURN(retval);
    
    /* Free all memory and data structures, EventSet must be empty. */
    if ( (retval=PAPI_destroy_eventset(&EventSet)) != PAPI_OK)
        ERROR_RETURN(retval);
    
    /* Free the resources used by PAPI */
    PAPI_shutdown();

    /* Free matrices */
    free_matrix(A,size);
    free_matrix(B,size);
    free_matrix(C,size);
    
    exit(0);
}

int atoloop_order(char* arg)
{
	if (strcmp(arg, "ijk") == 0)
	{
		return IJK;
	}
	else if (strcmp(arg, "jik") == 0)
	{
		return JIK;
	}
	else if (strcmp(arg, "jki") == 0)
	{
		return JKI;
	}
	else if (strcmp(arg, "kji") == 0)
	{
		return KJI;
	}
	else if (strcmp(arg, "ikj") == 0)
	{
		return IKJ;
	}
	else if (strcmp(arg, "kij") == 0)
	{
		return KIJ;
	}
	
	return ERR;
}

double** alloc_matrix(int size)
{
    int i;
    double** M = (double**) malloc(size * sizeof(double *));
    
    for (i=0; i<size; i++)
    {
        M[i] = (double*) malloc(size * sizeof(double));
    }
    
    return M;
}

void cache_reset()
{
	double values[32000];

	int i;
	
	for(i = 0; i < 32000; i++)
	{
		values[i] = -40000.0;
	}
}

void free_matrix(double** M, int size)
{
    int i;
    
    for (i=0; i<size; i++) 
    {
        free(M[i]);
    }
    
    free(M);
}

void matrix_ones_init(double** M, int size)
{
    int i, j;
    
    for (i=0; i<size; i++)
    {
        for (j=0; j<size; j++)
        {
            M[i][j] = 1.0;
        }
    }
}

void matrix_rand_init(double** M, int size)
{
    int i, j;
    
    for (i=0; i<size; i++)
    {
        for (j=0; j<size; j++)
        {
            M[i][j] = rand()/123.0;
        }
    }
}

void matrix_zero_init(double** M, int size)
{
    int i, j;
    
    for (i=0; i<size; i++)
    {
        for (j=0; j<size; j++)
        {
            M[i][j] = 0.0;
        }
    }
}

void print_matrix(double** M, int size)
{
    int i, j;
    
    for (i=0; i<size; i++)
    {
        for (j=0; j<size; j++)
        {
            printf("%f\t", M[i][j]);
        }

        printf("\n");
    }
    
    printf("\n");
}
